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ABSTRACT 



The Fisher information matrix of the cosmic microwave background (CMB) radiation power spectrum coefficients is a 
fundamental quantity that specifies the information content of a CMB experiment. In the most general case, its exact 
calculation scales with the third power of the number of data points A'' and is therefore computationally prohibitive for 
state-of-the-art surveys. Applicable to a very large class of CMB experiments without special symmetries, we show how 
to compute the Fisher matrix in only C(A'^^ log A'^) operations as long as the inverse noise covariance matrix can be 
applied to a data vector in time ©(^^^x log^max) . This assumption is true to a good approximation for all CMB data 
sets taken so far. The method takes into account common systematics such as arbitrary sky coverage and realistic noise 
correlations. As a consequence, optimal quadratic power spectrum estimation also becomes feasible in C'(A'^logA) 
operations for this large group of experiments. We discuss the relevance of our findings to other areas of cosmology 
where optimal power spectrum estimation plays a role. 
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1. Introduction 

Observations of the cosmic microwave background (CMB) 
radiation have proven to be a cornerstone of mod- 
ern cosmology, leading to a vigorou s experimental ef- 
fort to measure its anisotropies (e.g., Smoot et al. '199^ 
Netterfieldet al. 2002; Bennett et al. 2003; Kuo et al. 20o1 
bunklev et 811120 1 It iKeisler et al.ll2011l) . Whilst the limited 



texts, e.g., galaxy surveys (e.g.. lHamilton &: Tegmar3l2000l : 
|Seo fc Eiscnstein . ,2003 | ) . grav itational wave astronomy (e.g 



Be rti et all 120051: IVallisneril I2008D. weak l ensing surveys 
(e.g.. iHu fc TegmarkI llQOalKitching et al.' '200 3), and to 



angular resolution of the first CMB experiments naturally 
restricted the data size, observations obtained with Planck, 
the third generation CMB satellite experiment, and high- 
resolution ground-based experi ments will soon deliver maps 
with as many as 0(10'') pixels (jPlanck Collaboration et akl 
1201 ll) . challenging the performance of data analysis tools. 

The information on the most fundamental cosmologi- 
cal parameters given by a CMB sky map is fully contained 
in the much smaller set of angular power spe ctrum coeffi- 
cients (jJungman et al.lll996HBond et al.lll997D . This makes 
the power spectrum a convenient intermediate stage prod- 
uct in the analysis chain. For the lossless extraction of the 
parameters in a subsequent step, however, a thorough char- 
acterization of the statistical properties of the power spec- 
trum coefficients is necessary. 

As one of the most important objects in statistics, the 
Fisher information matrix quantifies the ability to constrain 
a set of parameters by means of an experiment. In addition, 
it reflects the (possibly complicated) correlation structure 
among them. Though formally only applicable to (multi- 
variate) Gaussian variables, the Fisher information matrix 
has been put to use in a wide variety of cosmological con- 



optimize numerical qutd^TtAre {e.g^Emith & Zaidarriagal 
120 Uh . etc. 

Unfortunately, the exact calculation of the Fisher ma- 
trix for the CMB power spectrum of present-day experi- 
ments is numerically intractable: in the most general case, 
the time complexity of suitable algorithms is O(iV^) , where 
N is the number of pixels of the CMB data map (jBorrilll 
IT999I) . To overcome this problem, approximate methods 
have been proposed to speed up the calculation. These 
methods may rely on the use of Monte Carl o averages to 
estimate expectation values ()Tegmar"3 11997), or numeri- 
cal d ifferentiation of the likelihood function ( Perotto et al.l 
I2OO6D . However, they are plagued by stability issues and 
unable to calculate off-diagonal elements to reasonable pre- 
cision. Thus, an exact scheme to evaluate the Fisher matrix 
with moderate computational cost has remained elusive un- 
til now. 

The paper is organized as follows. In Sect. [H we pro- 
vide a short overview of the mathematical framework that 
our approach is based on. In Sect. |3l we then introduce an 
efficient way of calculating the Fisher information matrix. 
We first consider experiments with a simplified scanning 
strategy to introduce the main ideas, and then generalize 
our findings to more generic CMB experiments. Finally, we 
summarize our results in Sect. HI 
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2. The likelihood function on the ring torus 

For an arbitrary CMB experiment with Gaussian noise, the 
hkehhood function £ of the data d is given by 



C{Ct\d) 



1 



V|2^(S(Q)+N)| 



X exp 



-\s{s{Ct)+^r^d 



(1) 



where we have introduced the noise covariance matrix, N, 
and the signal covariance matrix, S, a function of the CMB 
power spectrum coefficients Ci. 

For an efficient yet still exact evaluation of the likeli- 
hood, we start out by considering a CMB experiment that 
scans the sky on iso-latitude circles. The foundations of 
this m ethod are explained in detail in IWandelt fc Hans"enl 
(|2003( ). Imposing this restriction on the scanning strategy 
allows us to map the resulting time-ordered data (TOD) 
structure onto the ring torus. Owing to the periodic bound- 
ary conditions, this method enables us to work in the 
Fourier basis, where the algebraic expressions are simpler. 

Quantifying the signal covariance structure in Fourier 
space, we first write the noise-free sky temperature in terms 
of spherical harmonic coefficients of the signal map, a^m. 



\Xp. 



(2) 



where the index p runs over the Fourier modes in the direc- 
tion of the rings, and r specifies the index for the cross-ring 
direction. In Eq. we apply the definition of the Wigner 
rotation matrix to introduce the real quantity d 



,,0,0i) = e-'™^Ml„,(0)e-™'*S (3) 



where we choose 9 — 9s, the constant latitude of the exper- 
iment's spin axis. We also make use of the rotated beam X 
according to 



Xf 



2i- 



4tt 



(4) 



where the bim are the expansion coefficients of the beam 
function in spherical harmonics, and the Wigner small d 
matrix is evaluated at the opening angle of the scanning 
circles, 9o. We note that this framework allows an exact 
treatment of arbitrarily shaped beam patterns. The signal 
covariance matrix {T^^T^^,*,) is then given by 



drr'N 



dj.,p,Xl 



(5) 



where the Kronecker delta ensures the block diagonal struc- 
ture of the signal covariance matrix in Fourier space. 

After having specified the signal properties, we now 
quantify the noise correlations. Assuming a stationary 
Gaussian process, we can describe the noise properties in 
the TOD domain by a power spectrum P{k). We stress that 
this ansatz enables the exact treatment of correlated noise, 
a common systematic effect in real- world experiments. The 



noise covariance matrix {T^^T^*,) simplifies to 



e "P 

P m^m'—O 

X J2 e^^'^'^C(A,TO,m'), (6) 



A=0 



where is the number of rings in the data set, and A^p the 
number of pixels per ring. We have introduced an auxiliary 
function C(A,m,m'), which is defined as 



C(A,m,m') = P{k) 



e « 



fe(Af,.A-|-m-?Ti') 



(7) 



fc=0 



In analogy to S, we also find the noise covariance matrix 
N to be block diagonal in Fourier space. 

3. Calculating the Fisher matrix on the ring torus 

We now propose a strategy for the exact calculation of the 
Fisher matrix in only ©(A^^logA^) operations. Given the 
definition of the Fisher matrix as the covariance of the score 
function, from Eq. ([T]), we obtain (Tcg mark.1997.) 

52 \nC 



dCf , dCf . 



= -trfC-ip^^C-ip^ 
2 L 



(8) 



where C = S N, and = dC/dCi. The evaluation 
of Eq. ^ in its genera l form written here takes 0[N^) 
operations (|Borrilllll999l) . 

3.1. Idealized ring-torus scanning strategy 

We first consider an idealized experiment as described in 
the last section. For the calculation of the Fisher matrix, 
we need to compute the inverse of the covariance matrix, 
C~^, which can be done brute force in ^^(A^^) operations 
owing to its block diagonal structure. The second ingredi- 
ent, according to Eq. ([5]), is the derivative of the covariance 
matrix with respect to the power spectrum coefficients. On 
the ring torus, we can calculate the derivative for each block 
r separately as 



dCr 

da 



pp' 



dCe 



pp' 



— N'^ d^pXip dlpiXgp, 

Hrp Hrp' ' 



(9) 



that is, each block of the matrix is a rank one object and 
can therefore be written as the outer product of a single vec- 
tor q. Introducing the decomposition P^ — q^q^^, Eq. ^ 
simplifies considerably to 

r 

= lY.\^i^'c;\i^\'. (10) 
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Fig. 1. Fast calculation of the Fisher matrix. For several 
values of ^max, we show the measured runtime of our algo- 
rithm {filled circles). To guide the eye, we overplot a ^^ax 
power- law {dashed line). 

This exact rank-1 representation of the derivative of the sig- 
nal covariance matrix in the ring-torus Fourier domain is 
the key ingredient to our fast algorithm. We note that the 
analogous derivative in the spherical harmonic domain is 
also b lock diagonal but with blocks of rank 2i+l (jTegmarkl 
[l993)- We now show how to exploit this low-rank represen- 
tation. 

For each of the O(^max) independent blocks of the co- 
variance matrix, we perform the following calculations. We 
first evaluate the -^max matrix-vector products C~^q, at a 
numerical cost of C(^max)- We then accumulate each of 
the ^^ax entries of the Fisher matrix by a simple vector 
dot product, with a ©(^max) scaling. The overall time com- 
plexity of the algorithm therefore amounts to C(^max) ~ 
O(iV^) operations. 

To confirm the scaling behavior as claimed above, we 
did a full numerical implementation of the algorithm. In 
Fig.[Tl we show the timings we obtained for different values 
of ^max- The measured runtime clearly follows the predicted 
time complexity O(^max)- The tests were carried out on a 
cluster of Intel Xenon processors with a 3 GHz clock rate 
using 256 CPU cores. 

3.2. Generalization beyond ring-torus scans and arbitrary sky 
coverage 

After having outlined the algorithm for the specific case of 
a simplified scanning strategy, we now discuss how to gener- 
alize the results to a generic CMB experiment. To this end, 
we propose a resampling of the input data to a grid that can 
be mapped onto the ring torus (e.g., an ECP grid, where 
Os = do = 7r/2). This has the advantage that the simpli- 
fying signal covariance properties, as described above, are 
still fulfilled, i.e. Eq. 1^ still holds. Unfortunately, as the 



grid no longer follows the scanning strategy, the possibility 
of achieving an exact treatment of asymmetric beams at 
no additional computational costs is lost. In addition, the 
noise properties may be more complicated, owing to, e.g., 
anisotropically distributed observation time, or the pres- 
ence of a Galactic mask. As a result, an explicit expression 
for the covariance matrix C may not exist. Luckily, we only 
have to evaluate products of inverse covariance matrix and 
sky map. This problem is standard in CMB data analysis, 
and iterative solvers have been successfully developed for 
this problem (jSmith et al.ll2007[ ). 

As the block diagonal structure of C is destroyed, we 
solve 

FeA = l\qlC-\eJ\ (11) 

by first finding solutions to the ^max equations x — C^^q 
iteratively, and then contracting the outcome with a sim- 
ple dot product. Once the calculation has been completed 
successfully, the resulting Fisher matrix is independent of 
the pixelization used for the actual computation. 

The overall 0{N^) time complexity of the algorithm 
can be sustained, if the solver scales as 0{£^^^). In the 
most general case, where the full characterization of the 
noise properties requires a dense N x N noise covariance 
matrix, the numerical costs increase to C(£^ax) ~ 0[N^^^) 
operations. Even worse, the memory requirement to store 
the matrix rises to 0{N'^), reaching thousands of terabytes 
for Planck-sized data. 

Luckily, all carefully designed CMB experiments allow 
either a structured or a sparse representation of the noise 
covariance matrix. In the presence of an anisotropic yet 
uncorrelated noise and arbitrary sky coverage, for exam- 
ple, it reduces to a diagonal matrix in a real space rep- 
resentation. This ansatz made an ef ficient inverse varianc e 
weighting of WMAP data possible (jKomatsu et al.ll201lh . 
Under very general circumstances, the noise covariance can 
be written in a mixed real-space and time-ordered domain 
representation, e.g. N^^ = A^N^q^A, where A is the 
pointing matrix and Ntod is the noise covariance in the 
time-ordered domain. For essentially all past and current 
CMB experiments, the noise is piecewise stationary in the 
time-ordered domain, which means that Ntod is block- 
Toeplitz, and the matrix-vector product of with the 
data takes 0(^^x00 log A^tod) operations. Since there are 
at most C'(^max) orientations in which a CMB experiment 
can scan the sky, A^tod ^ ^max- ^s a consequence, the cal- 
culation of the Fisher matrix in 0{N^ log N) operations 
becomes feasible. 

A straightforward treatment of asymmetric beams is 
only possible if the data grid follows the scanning strategy 
directly. Although this is not the case in the general setup 
discussed above, we note that asymmetric beams can be in- 
cluded into the analysis at additional computational costs 
of C(™max)j determined by the azimuthal structure of the 
beams. For example, the treatment of elliptical Gaussian 
beams increases the numerical complexity by a factor of 
four. 

4. Discussion and conclusions 

The Fisher information matrix is a fundamental quantity 
in statistics, reflecting the predictive power of the experi- 
ment under consideration. Despite its importance, thorough 



3 



Franz Eisner and Benjamin D. Wandelt: Fast calculation of the Fisher matrix 



systematic Fisher-matrix-based studies of experiments, de- 
signed to measure the CMB power spectrum, have not yet 
been conducted. The reason for this shortcoming can be 
found in the associated computational expenses: general 
algorithms show a time complexity of 0(^N^), where N is 
the number of pixels of the survey, renderi ng an analy sis 
for state of the art experiments impossible (|Borrilllll999[ ). 

Here, we have presented a new method for the exact cal- 
culation of the Fisher matrix of the CMB power spectrum 
coefficients in only O^N^logN) operations. To do so, we 
first considered experiments with a specific scanning strat- 
egy, where the sky is observed on iso-latitude circles. This 
restriction has allowed us to map the TOD onto the ring 
torus in order to take advantage of the symmetries of this 
manifold. We then cast the equation for the CMB likeli- 
hood function into Fourier space and found that both the 
signal and the noise covar iance matrices are block diagonal 
(jWandelt &: Hansenll2003f) . For each block, the derivative of 
the signal covariance matrix with respect to a single power 
spectrum coefhcient was also found to be a rank one ma- 
trix. As a result, we made significant progress in ensuring 
that the equations to calculate the Fisher matrix exactly 
simplify, requiring only ©(iV^) operations to evaluate. 

We then relaxed the assumption of a simplified scan- 
ning strategy. If iterative methods allowed us to calculate 
inverse- variance- weighted sky maps in 0(^max log^max) op- 
erations, we showed that a 0[N^ log N) time complexity of 
our algorithm could be realized. To this end, the data were 
resampled onto an ECP grid, which is homeomorphic to the 
ring torus, where the properties of the signal covariance ma- 
trix were still simplified. Our fast method is applicable to 
a very general class of experiments, including, to a good 
approximation, all CMB data sets taken up to now. 

We consider a fast means of calculating the Fisher infor- 
mation matrix for CMB experiments of great importance. 
This new method not only allows us to study quantita- 
tively the impact of the common systematics of realistic 
CMB experiments, such as partial sky coverage, asymmet- 
ric beams, or correlated noise, on the power spectrum. The 
Fisher matrix also appears as a normalization fact or for un- 
biase d and optimal power spectrum estimators ()Tegmarkl 
|1997( ). Since its fast calculation has not been possible, until 
now, only approximative pseudo-Cf estimators have been 
efficient enough to be applicable to more recent CMB ex- 
periments (e.g.. ISzapudi et al.l 120011: IWandelt et"all 120011: 
iHivon et al.| l2002r The new scheme outlined here has the 
potential to lift this restriction. 

For the presentation of our method, we explicitly con- 
sidered an experiment designed to measure the CMB power 
spectrum. However, the described framework is general 
enough to be of value in other fields of application, where 
the statistical properties of an isotropic signal in spheri- 
cal geometry are investigated. These applications may in- 
clude, but are not limited to, the analysis of gala xy an- 
gular power spectra (Rim es fc Hamilton ,2005.: Lee fc~Penl 
2008^■ or weak lensiii g tomography ( Bernstein fc JainI 2004 
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